clear all
capture log close
capture set virtual off
set memory 1000m
capture set matsize 100000000000
capture set more off
capture virtual off
memory

set type double
set more off
set maxvar 32767
pause on
version 12.1

* <><><><><><><><><> *
* Mandy's computer 
* <><><><><><><><><> *

cd "/Users/Mandywu/Dropbox (MIT)/PNAS revision/code"


********** 
use campaign_all_10wk_nt, replace
bysort campaign_treated campaigntime: egen avglog_concentration=mean(log_concentration_SO2)
gen log_concentration_control = avglog_concentration if campaign_treated == 0
gen log_concentration_treated = avglog_concentration if campaign_treated == 1

drop if campaigntime == 15

rename log_concentration_treated treated_cities
rename log_concentration_control control_cities

generate age1 = (campaigntime + 5) 
replace  age1 = 0 if campaigntime > -5
generate age2 = (campaigntime + 5) 
replace  age2 = 0 if campaigntime < -4 | campaigntime > -1
generate age3 = (campaigntime + 5) 
replace  age3 = 0 if campaigntime < 0 | campaigntime2 > 0
generate age4 = (campaigntime + 5) 
replace  age4 = 0 if campaigntime2 < 1

generate int1 = 0
replace  int1 = 1 if campaigntime <= -5
generate int2 = 0
replace  int2 = 1 if campaigntime >= -4 & campaigntime <= -1
generate int3 = 0
replace  int3 = 1 if campaigntime >= 0 & campaigntime2 <= 0
generate int4 = 0
replace  int4 = 1 if campaigntime2 >= 1

regress treated_cities int1 int2 int3 int4 age1 age2 age3 age4, hascons
predict yhat_treated

regress control_cities int1 int2 int3 int4 age1 age2 age3 age4, hascons
predict yhat_control

	 
twoway (scatter treated_cities control_cities campaigntime, msize(3pt 3pt)) ///
         (line yhat_treated yhat_control campaigntime if campaigntime <= -5, sort lpattern(solid shortdash) lwidth(medium medthick) lcolor(navy maroon)) ///
		 (line yhat_treated yhat_control campaigntime if campaigntime >= -4 & campaigntime <= -1, sort lpattern(solid shortdash) lwidth(medium medthick) lcolor(navy maroon)) ///
		 (line yhat_treated yhat_control campaigntime if campaigntime >= 0 & campaigntime2 <= 0, sort lpattern(solid shortdash) lwidth(medium medthick) lcolor(navy maroon)) ///
		 (line yhat_treated yhat_control campaigntime if campaigntime2 >= 1, sort lpattern(solid shortdash) lwidth(medium medthick) lcolor(navy maroon)), ///
		 xline(-5, lpattern(solid) lcolor(navy)) ///
		 graphregion(color(white)) bgcolor(white) xlabel(-10 "-10" -9 " " -8 " " -7 " " -6 " " -5 " " -4 " " ///
		 -3 " " -2 "Annouce" -1 " " 0 "0" 1 " " 2"On-site" 3 " " 4 " " 5 " " 6 "+1" 7 " " 8 " " 9 " " 10 "+5" 11 " " 12 " " 13 " " 14 " " 15 "+10") ///
		  xline(-4.5 -0.5, lcolor(gs15)) ///
		 xline(-4, lwidth(6.5) lc(gs15)) xline(-3, lwidth(6.5) lc(gs15)) xline(-2, lwidth(6.5) lc(gs15)) xline(-1, lwidth(6.5) lc(gs15)) ///
		 xline(1 5, lcolor(gs12)) xline(1.5, lwidth(6.5) lc(gs12)) xline(2.5, lwidth(6.5) lc(gs12)) xline(3.5, lwidth(6.5) lc(gs12)) ///
		 xline(4.5, lwidth(6.5) lc(gs12)) xline(0, lpattern(solid) lc(edkblue)) ///
		 xtitle("Week relative to central inspection") ///
		 xlabel(,labsize(medium)) ytitle("log(SO2)", size(3)) graphregion(fcolor(white)) legend(order(1 "Plants in treated cities" 2 "Plants in control cities"))
graph export figureS2.png, replace



***************** figure S3
*****************
use campaign_all_10wk_match_nt, replace

loc minbalancedeventtime = -10
loc maxbalancedeventtime = 10

forval v = `minbalancedeventtime'/-1 {
loc absv = abs(`v')
*g d_event_neg`absv' = eventtime==`v' 
g d_event_neg`absv'_tr = campaigntime==`v' & campaign_treated == 1
}

g d_event_exact`0'_1tr = campaigntime==0 & campaign_treated == 1
g d_event_exact`0'_2tr = campaigntime==1 & campaign_treated == 1
g d_event_exact`0'_3tr = campaigntime==2 & campaign_treated == 1
g d_event_exact`0'_4tr = campaigntime==3 & campaign_treated == 1
g d_event_exact`0'_5tr = campaigntime==4 & campaign_treated == 1
g d_event_exact`0'_6tr = campaigntime==5 & campaign_treated == 1 if province_name == "河南省" | province_name == "黑龙江省" | province_name == "广东省" | province_name == "湖北省"

forval v = 1/`maxbalancedeventtime' {
*g d_event_pos`v' = eventtime==`v'
g d_event_pos`v'_tr = campaigntime2==`v' & campaign_treated == 1
}

***** add in all the result
sort plant_name yweek
by plant_name: egen weekmin= min(yweek)
gen t= yweek - weekmin
gen t2=t^2
save event_10wk, replace

**** anohter set of result - only change to -12/12

forval o=0/2 {
eststo y`o':  reghdfe log_concentration_SO2 d_event_neg10_tr-d_event_neg6_tr d_event_neg4_tr-d_event_pos10_tr [pweight=_webal], absorb(plant_name1 yweek) vce(cluster city_name)
}
eststo compaign_10wk_city


coefplot compaign_10wk_city, drop(campaign yweek *.month *.year _cons other_period) ///
graphregion(color(white)) bgcolor(white) xlabel(1 "-10" 2 " " 3 " " 4 " " 5 " " 6 " " 7 " " ///
8 "annouce" 9 " " 9.5 "0" 10 " " 11 " " 12 " " 13 "onsite" 14 " " 15 " " 16 "+1" 17 " " 18 " " 19 "+4" 20 " " 21 " " ///
22 "+7" 23 " " 24 " " 25 "+10" ) xline(9.5, lpattern(solid) lc(dkgreen)) xline(5.5, lpattern(dash)) ///
yline(0, lcolor(black)) xaxis(1 2) xline(6 9, lcolor(gs15)) xline(6.5, lwidth(6.5) lc(gs15)) xline(7.5, lwidth(6.5) lc(gs15)) ///
xline(8.5, lwidth(6.5) lc(gs15)) xline(10 15, lcolor(gs12)) xline(10.5, lwidth(6.5) lc(gs12)) ///
xline(11.5, lwidth(6.5) lc(gs12)) xline(12.5, lwidth(6.5) lc(gs12)) xline(13.5, lwidth(6.5) lc(gs12)) ///
xline(14.5, lwidth(6.5) lc(gs12)) title("Centered around -10/+10 weeks of arrival") /// 
xtitle("Week relative to onsite inspection", axis(1)) xtitle("", axis(2)) ///
vertical msymbol(s) ciopts(recast(rcap)) xlabel(,labsize(small)) ytitle("log(SO2)", size(3)) graphregion(fcolor(white)) baselevels
graph export figureS3_nt.png, replace



*****************
use campaign_all_10wk_match, replace

loc minbalancedeventtime = -10
loc maxbalancedeventtime = 10

forval v = `minbalancedeventtime'/-1 {
loc absv = abs(`v')
*g d_event_neg`absv' = eventtime==`v' 
g d_event_neg`absv'_tr = campaigntime==`v' & campaign_treated == 1
}

g d_event_exact`0'_1tr = campaigntime==0 & campaign_treated == 1
g d_event_exact`0'_2tr = campaigntime==1 & campaign_treated == 1
g d_event_exact`0'_3tr = campaigntime==2 & campaign_treated == 1
g d_event_exact`0'_4tr = campaigntime==3 & campaign_treated == 1
g d_event_exact`0'_5tr = campaigntime==4 & campaign_treated == 1
g d_event_exact`0'_6tr = campaigntime==5 & campaign_treated == 1 if province_name == "河南省" | province_name == "黑龙江省" | province_name == "广东省" | province_name == "湖北省"

forval v = 1/`maxbalancedeventtime' {
*g d_event_pos`v' = eventtime==`v'
g d_event_pos`v'_tr = campaigntime2==`v' & campaign_treated == 1
}

***** add in all the result
sort plant_name yweek
by plant_name: egen weekmin= min(yweek)
gen t= yweek - weekmin
gen t2=t^2
save event_10wk, replace

**** anohter set of result - only change to -12/12

forval o=0/2 {
eststo y`o':  reghdfe log_concentration_SO2 d_event_neg10_tr-d_event_neg6_tr d_event_neg4_tr-d_event_pos10_tr [pweight=_webal], absorb(plant_name1 yweek) vce(cluster city_name)
}
eststo compaign_10wk_city


coefplot compaign_10wk_city, drop(campaign yweek *.month *.year _cons other_period) ///
graphregion(color(white)) bgcolor(white) xlabel(1 "-10" 2 " " 3 " " 4 " " 5 " " 6 " " 7 " " ///
8 "annouce" 9 " " 9.5 "0" 10 " " 11 " " 12 " " 13 "onsite" 14 " " 15 " " 16 "+1" 17 " " 18 " " 19 "+4" 20 " " 21 " " ///
22 "+7" 23 " " 24 " " 25 "+10" ) xline(9.5, lpattern(solid) lc(dkgreen)) xline(5.5, lpattern(dash)) ///
yline(0, lcolor(black)) xaxis(1 2) xline(6 9, lcolor(gs15)) xline(6.5, lwidth(6.5) lc(gs15)) xline(7.5, lwidth(6.5) lc(gs15)) ///
xline(8.5, lwidth(6.5) lc(gs15)) xline(10 15, lcolor(gs12)) xline(10.5, lwidth(6.5) lc(gs12)) ///
xline(11.5, lwidth(6.5) lc(gs12)) xline(12.5, lwidth(6.5) lc(gs12)) xline(13.5, lwidth(6.5) lc(gs12)) ///
xline(14.5, lwidth(6.5) lc(gs12)) title("Centered around -10/+10 weeks of arrival") /// 
xtitle("Week relative to onsite inspection", axis(1)) xtitle("", axis(2)) ///
vertical msymbol(s) ciopts(recast(rcap)) xlabel(,labsize(small)) ytitle("log(SO2)", size(3)) graphregion(fcolor(white)) baselevels
graph export figureS3_all.png, replace


*****************
use campaign_all_10wk_match_at, replace

loc minbalancedeventtime = -10
loc maxbalancedeventtime = 10

forval v = `minbalancedeventtime'/-1 {
loc absv = abs(`v')
*g d_event_neg`absv' = eventtime==`v' 
g d_event_neg`absv'_tr = campaigntime==`v' & campaign_treated == 1
}

g d_event_exact`0'_1tr = campaigntime==0 & campaign_treated == 1
g d_event_exact`0'_2tr = campaigntime==1 & campaign_treated == 1
g d_event_exact`0'_3tr = campaigntime==2 & campaign_treated == 1
g d_event_exact`0'_4tr = campaigntime==3 & campaign_treated == 1
g d_event_exact`0'_5tr = campaigntime==4 & campaign_treated == 1
g d_event_exact`0'_6tr = campaigntime==5 & campaign_treated == 1 if province_name == "河南省" | province_name == "黑龙江省" | province_name == "广东省" | province_name == "湖北省"

forval v = 1/`maxbalancedeventtime' {
*g d_event_pos`v' = eventtime==`v'
g d_event_pos`v'_tr = campaigntime2==`v' & campaign_treated == 1
}

***** add in all the result
sort plant_name yweek
by plant_name: egen weekmin= min(yweek)
gen t= yweek - weekmin
gen t2=t^2
save event_10wk, replace

**** anohter set of result - only change to -12/12

forval o=0/2 {
eststo y`o':  reghdfe log_concentration_SO2 d_event_neg10_tr-d_event_neg6_tr d_event_neg4_tr-d_event_pos10_tr [pweight=_webal], absorb(plant_name1 yweek) vce(cluster city_name)
}
eststo compaign_10wk_city


coefplot compaign_10wk_city, drop(campaign yweek *.month *.year _cons other_period) ///
graphregion(color(white)) bgcolor(white) xlabel(1 "-10" 2 " " 3 " " 4 " " 5 " " 6 " " 7 " " ///
8 "annouce" 9 " " 9.5 "0" 10 " " 11 " " 12 " " 13 "onsite" 14 " " 15 " " 16 "+1" 17 " " 18 " " 19 "+4" 20 " " 21 " " ///
22 "+7" 23 " " 24 " " 25 "+10" ) xline(9.5, lpattern(solid) lc(dkgreen)) xline(5.5, lpattern(dash)) ///
yline(0, lcolor(black)) xaxis(1 2) xline(6 9, lcolor(gs15)) xline(6.5, lwidth(6.5) lc(gs15)) xline(7.5, lwidth(6.5) lc(gs15)) ///
xline(8.5, lwidth(6.5) lc(gs15)) xline(10 15, lcolor(gs12)) xline(10.5, lwidth(6.5) lc(gs12)) ///
xline(11.5, lwidth(6.5) lc(gs12)) xline(12.5, lwidth(6.5) lc(gs12)) xline(13.5, lwidth(6.5) lc(gs12)) ///
xline(14.5, lwidth(6.5) lc(gs12)) title("Centered around -10/+10 weeks of arrival") /// 
xtitle("Week relative to onsite inspection", axis(1)) xtitle("", axis(2)) ///
vertical msymbol(s) ciopts(recast(rcap)) xlabel(,labsize(small)) ytitle("log(SO2)", size(3)) graphregion(fcolor(white)) baselevels
graph export figureS3_at.png, replace


****** figure S4
use panel_plant_weekly_diff_control, replace

encode plant_name, generate(plant_name1)
encode city_name, generate(city_name1)
gen campaign_date = 0

replace campaign_date = yweek if (first_period==1 & (province_name == "河南省" | province_name == "黑龙江省")) | ///
(second_period==1 & (province_name == "广东省" | province_name=="湖北省")) | ///
(third_period==1 & (province_name == "山西省" | province_name =="湖南省")) | ///
(fourth_period==1  & (province_name == "山东省" | province_name=="浙江省"))

sort plant_name yweek
by plant_name: gen weeknum1=_n
by plant_name: gen target1=weeknum1 if yweek==campaign_date
egen td1=min(target1), by(plant_name)
gen campaigntime=weeknum1-td1
egen td2=max(target1), by(plant_name)
gen campaigntime2=weeknum1-td2
*drop weeknum1 target1 td1 td2
gen baseline = td1 - 5
rename td1 first_treat

*keep if campaigntime >= -30 & campaigntime2 <= 30
keep if year >= 2016 & year <= 2017

save campaign_est_baseline.dta, replace

*** here the treatment week is the first week 
csdid  log_concentration_SO2 per_capita_log, ivar(plant_name1) time(yweek) gvar(baseline) method(dripw) cluster(city_name1)
*estat event, window(-10 14) estore(eq1)
estat event, window(-5 19) estore(eq1)
esttab eq1 using event_csdid_nocontrol.tex, se replace

event_plot eq1, default_look ciplottype(rcap)  alpha(0.1) ///
graph_opt(xtitle("Periods since the announcement") ytitle("Average causal effect") ///
title("Event study") xlabel(-5 "-10" -4 " " -3 " " -2 " " -1 "-6" 0 " " 1 " " 2 "Annouce" 3 " " 4 " " 5 " " 6 " " 7 "On-site" 8 " " 9 " " 10 "+1" 11 " " 12 " " 13 " " 14 "+5" 15 " " 16 " " 17 " " 18 " " 19 "+10") ///
xline(0 4, lcolor(gs15)) xline(0.5, lwidth(6.5) lc(gs15)) xline(1.5, lwidth(6.5) lc(gs15)) xline(2.5, lwidth(6.5) lc(gs15)) ///
xline(3.5, lwidth(6.5) lc(gs15)) xline(5 9, lcolor(gs12)) xline(5.5, lwidth(6.5) lc(gs12)) xline(6.5, lwidth(6.5) lc(gs12)) ///
xline(7.5, lwidth(6.5) lc(gs12)) xline(8.5, lwidth(6.5) lc(gs12)) ///
xline(4.5, lpattern(solid) lc(dkgreen)) xline(-0.5, lpattern(dash) lc(dkgreen)) yline(0, lpattern(dash) lcolor(black)) ///
xlabel(,labsize(small)) graphregion(fcolor(white))) stub_lag(Tp#) stub_lead(Tm#) together
graph export figureS4.png, replace


/***************** figure 4a
use henan_monitor, replace
encode plant_name, generate(plant_name1)
keep if province_name == "河南省"

gen event_date = 0
replace event_date = yweek if inspection==1
sort plant_name yweek
by plant_name: gen weeknum=_n
by plant_name: gen target=weeknum if yweek==event_date
egen td=min(target), by(plant_name)
drop target
gen eventtime=weeknum-td

keep if eventtime>=-12 & eventtime <= 12
bysort treated eventtime: egen avglog_concentration_SO2=mean(log_concentration_SO2)
gen log_concentration_SO2_control = avglog_concentration_SO2 if treated == 0
gen log_concentration_SO2_treated = avglog_concentration_SO2 if treated == 1

rename log_concentration_SO2_treated treated_power_plants
rename log_concentration_SO2_control control_power_plants

line treated_power_plants control_power_plants eventtime, ///
graphregion(color(white)) bgcolor(white) xlabel(-12 "-12" -11 " " -10 " " -9 "-9" -8 " " -7 " " -6 "prepare" -5 " " -4 " " ///
-3 "-3" -2 " " -1 " " 0 "onsite" 1 " " 2" " 3 "3" 4 " " 5 " " 6 "6" 7 " " 8 " " 9 "9" 10 " " 11 " " 12 "12" ) ///
 yline(0, lcolor(black)) xaxis(1 2) xline(-8, lwidth(5.5) lc(gs15)) xline(-7, lwidth(5.5) lc(gs15)) ///
xline(-6, lwidth(5.5) lc(gs15)) xline(-5, lwidth(5.5) lc(gs15)) xline(-4, lwidth(5.5) lc(gs15)) xline(-3, lwidth(5.5) lc(gs15)) ///
xline(-2.5 2.5, lcolor(gs12)) xline(-2, lwidth(5.5) lc(gs12)) xline(-1, lwidth(5.5) lc(gs12)) xline(0, lwidth(5.5) lc(gs12)) ///
xline(1, lwidth(5.5) lc(gs12)) xline(2, lwidth(5.5) lc(gs12)) xline(-9, lpattern(dash)) xline(0, lpattern(solid) lc(dkgreen)) ///
xla(-9 "tips" 0 "complaints", axis(2) labsize(medium)) xtitle("Week relative to being complained", axis(1)) xtitle("", axis(2)) ///
xlabel(,labsize(medium)) ytitle("log(SO2) ug/m^3", size(3)) graphregion(fcolor(white)) legend(order(1 "Treated power plants" 2 "Control power plants"))
graph export figureS3_a.png, replace

******************** figure 4b
use henan_cems, replace
encode plant_name, generate(plant_name1)

gen log_concentration_SO2 = log(concentration_SO2)

gen event_date = 0
replace event_date = yweek if inspection==1
sort plant_name yweek
by plant_name: gen weeknum=_n
by plant_name: gen target=weeknum if yweek==event_date
egen td=min(target), by(plant_name)
drop target
gen eventtime=weeknum-td

keep if eventtime>=-12 & eventtime <= 12
bysort treated eventtime: egen avglog_concentration_SO2=mean(log_concentration_SO2)
gen log_concentration_SO2_control = avglog_concentration_SO2 if treated == 0
gen log_concentration_SO2_treated = avglog_concentration_SO2 if treated == 1

rename log_concentration_SO2_treated treated_power_plants
rename log_concentration_SO2_control control_power_plants

line treated_power_plants control_power_plants eventtime, ///
graphregion(color(white)) bgcolor(white) xlabel(-12 "-12" -11 " " -10 " " -9 "-9" -8 " " -7 " " -6 "prepare" -5 " " -4 " " ///
-3 "-3" -2 " " -1 " " 0 "onsite" 1 " " 2" " 3 "3" 4 " " 5 " " 6 "6" 7 " " 8 " " 9 "9" 10 " " 11 " " 12 "12" ) ///
 yline(0, lcolor(black)) xaxis(1 2) xline(-8, lwidth(5.5) lc(gs15)) xline(-7, lwidth(5.5) lc(gs15)) ///
xline(-6, lwidth(5.5) lc(gs15)) xline(-5, lwidth(5.5) lc(gs15)) xline(-4, lwidth(5.5) lc(gs15)) xline(-3, lwidth(5.5) lc(gs15)) ///
xline(-2.5 2.5, lcolor(gs12)) xline(-2, lwidth(5.5) lc(gs12)) xline(-1, lwidth(5.5) lc(gs12)) xline(0, lwidth(5.5) lc(gs12)) ///
xline(1, lwidth(5.5) lc(gs12)) xline(2, lwidth(5.5) lc(gs12)) xline(-9, lpattern(dash)) xline(0, lpattern(solid) lc(dkgreen)) ///
xla(-9 "tips" 0 "complaints", axis(2) labsize(medium)) xtitle("Week relative to being complained", axis(1)) xtitle("", axis(2)) ///
xlabel(,labsize(medium)) ytitle("log(SO2) mg/m^3", size(3)) graphregion(fcolor(white)) legend(order(1 "Treated power plants" 2 "Control power plants"))
graph export figureS3_b.png, replace



